#import imageio.v2 as imageio
import imageio
import pydicom
from pydicom.data import get_testdata_file
import numpy as np
# load a single DICOM image from the scan volume
ct_file = get_testdata_file("CT_small.dcm")
mr_file = get_testdata_file("MR_small.dcm")
im_ct = imageio.imread(ct_file, progress=True)
im_mr = imageio.imread(mr_file, progress=True)
C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/3685182357.py:11: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly. im_ct = imageio.imread(ct_file, progress=True) C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/3685182357.py:12: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly. im_mr = imageio.imread(mr_file, progress=True)
Comprising : Patient Demographics and Acquisition Information
ds_ct = pydicom.dcmread(ct_file)
ds_ct
Dataset.file_meta ------------------------------- (0002, 0000) File Meta Information Group Length UL: 192 (0002, 0001) File Meta Information Version OB: b'\x00\x01' (0002, 0002) Media Storage SOP Class UID UI: CT Image Storage (0002, 0003) Media Storage SOP Instance UID UI: 1.3.6.1.4.1.5962.1.1.1.1.1.20040119072730.12322 (0002, 0010) Transfer Syntax UID UI: Explicit VR Little Endian (0002, 0012) Implementation Class UID UI: 1.3.6.1.4.1.5962.2 (0002, 0013) Implementation Version Name SH: 'DCTOOL100' (0002, 0016) Source Application Entity Title AE: 'CLUNIE1' ------------------------------------------------- (0008, 0005) Specific Character Set CS: 'ISO_IR 100' (0008, 0008) Image Type CS: ['ORIGINAL', 'PRIMARY', 'AXIAL'] (0008, 0012) Instance Creation Date DA: '20040119' (0008, 0013) Instance Creation Time TM: '072731' (0008, 0014) Instance Creator UID UI: 1.3.6.1.4.1.5962.3 (0008, 0016) SOP Class UID UI: CT Image Storage (0008, 0018) SOP Instance UID UI: 1.3.6.1.4.1.5962.1.1.1.1.1.20040119072730.12322 (0008, 0020) Study Date DA: '20040119' (0008, 0021) Series Date DA: '19970430' (0008, 0022) Acquisition Date DA: '19970430' (0008, 0023) Content Date DA: '19970430' (0008, 0030) Study Time TM: '072730' (0008, 0031) Series Time TM: '112749' (0008, 0032) Acquisition Time TM: '112936' (0008, 0033) Content Time TM: '113008' (0008, 0050) Accession Number SH: '' (0008, 0060) Modality CS: 'CT' (0008, 0070) Manufacturer LO: 'GE MEDICAL SYSTEMS' (0008, 0080) Institution Name LO: 'JFK IMAGING CENTER' (0008, 0090) Referring Physician's Name PN: '' (0008, 0201) Timezone Offset From UTC SH: '-0500' (0008, 1010) Station Name SH: 'CT01_OC0' (0008, 1030) Study Description LO: 'e+1' (0008, 1090) Manufacturer's Model Name LO: 'RHAPSODE' (0009, 0010) Private Creator LO: 'GEMS_IDEN_01' (0009, 1001) [Full fidelity] LO: 'GE_GENESIS_FF' (0009, 1002) [Suite id] SH: 'CT01' (0009, 1004) [Product id] SH: 'HiSpeed CT/i' (0009, 1027) [Image actual date] SL: 862399669 (0009, 1030) [Service id] SH: '' (0009, 1031) [Mobile location number] SH: '' (0009, 10e6) [Genesis Version - now] SH: '05' (0009, 10e7) [Exam Record checksum] UL: 973283917 (0009, 10e9) [Actual series data time stamp] SL: 862399669 (0010, 0010) Patient's Name PN: 'CompressedSamples^CT1' (0010, 0020) Patient ID LO: '1CT1' (0010, 0030) Patient's Birth Date DA: '' (0010, 0040) Patient's Sex CS: 'O' (0010, 1002) Other Patient IDs Sequence 2 item(s) ---- (0010, 0020) Patient ID LO: 'ABCD1234' (0010, 0022) Type of Patient ID CS: 'TEXT' --------- (0010, 0020) Patient ID LO: '1234ABCD' (0010, 0022) Type of Patient ID CS: 'TEXT' --------- (0010, 1010) Patient's Age AS: '000Y' (0010, 1030) Patient's Weight DS: '0.0' (0010, 21b0) Additional Patient History LT: '' (0011, 0010) Private Creator LO: 'GEMS_PATI_01' (0011, 1010) [Patient Status] SS: 0 (0018, 0010) Contrast/Bolus Agent LO: 'ISOVUE300/100' (0018, 0022) Scan Options CS: 'HELICAL MODE' (0018, 0050) Slice Thickness DS: '5.0' (0018, 0060) KVP DS: '120.0' (0018, 0088) Spacing Between Slices DS: '5.0' (0018, 0090) Data Collection Diameter DS: '480.0' (0018, 1020) Software Versions LO: '05' (0018, 1040) Contrast/Bolus Route LO: 'IV' (0018, 1100) Reconstruction Diameter DS: '338.6716' (0018, 1110) Distance Source to Detector DS: '1099.3100585938' (0018, 1111) Distance Source to Patient DS: '630.0' (0018, 1120) Gantry/Detector Tilt DS: '0.0' (0018, 1130) Table Height DS: '133.699997' (0018, 1150) Exposure Time IS: '1601' (0018, 1151) X-Ray Tube Current IS: '170' (0018, 1152) Exposure IS: '170' (0018, 1160) Filter Type SH: 'LARGE BOWTIE FIL' (0018, 1190) Focal Spot(s) DS: '0.7' (0018, 1210) Convolution Kernel SH: 'STANDARD' (0018, 5100) Patient Position CS: 'FFS' (0019, 0010) Private Creator LO: 'GEMS_ACQU_01' (0019, 1002) [Detector Channel] SL: 912 (0019, 1003) [Cell number at Theta] DS: '373.75' (0019, 1004) [Cell spacing] DS: '1.0166' (0019, 100f) [Horiz. Frame of ref.] DS: '955.799988' (0019, 1011) [Series contrast] SS: 2 (0019, 1013) [Start number for baseline] SS: 0 (0019, 1014) [End number for baseline] SS: 0 (0019, 1015) [Start number for enhanced scans] SS: 0 (0019, 1016) [End number for enhanced scans] SS: 0 (0019, 1017) [Series plane] SS: 2 (0019, 1018) [First scan ras] LO: 'S' (0019, 1019) [First scan location] DS: '7.79187' (0019, 101a) [Last scan ras] LO: 'I' (0019, 101b) [Last scan loc] DS: '-320.197968' (0019, 101e) [Display field of view] DS: '0.0' (0019, 1023) [Table Speed [mm/rotation]] DS: '5.0' (0019, 1024) [Mid Scan Time [sec]] DS: '17.784578' (0019, 1025) [Mid scan flag] SS: 1 (0019, 1026) [Tube Azimuth [degree]] SL: 0 (0019, 1027) [Rotation Speed [msec]] DS: '1.0' (0019, 102a) [x-ray On position] DS: '178.079926' (0019, 102b) [x-ray Off position] DS: '3994.299316' (0019, 102c) [Number of triggers] SL: 10431 (0019, 102e) [Angle of first view] DS: '-718.079956' (0019, 102f) [Trigger frequency] DS: '984.0' (0019, 1039) [SFOV Type] SS: 16 (0019, 1040) [Stat recon flag] SS: 0 (0019, 1041) [Compute type] SS: 1 (0019, 1042) [Segment Number] SS: 0 (0019, 1043) [Total Segments Required] SS: 0 (0019, 1044) [Interscan delay] DS: '1.0' (0019, 1047) [View compression factor] SS: 1 (0019, 104a) [Total no. of ref channels] SS: 6 (0019, 104b) [Data size for scan data] SL: 10431 (0019, 1052) [Recon post proc. Flag] SS: 1 (0019, 1057) [CT water number] SS: -95 (0019, 1058) [CT bone number] SS: 0 (0019, 105e) [Number of channels (1...512)] SL: 763 (0019, 105f) [Increment between channels] SL: 1 (0019, 1060) [Starting view] SL: 1969 (0019, 1061) [Number of views] SL: 1576 (0019, 1062) [Increment between views] SL: 1 (0019, 106a) [Dependent on #views processed] SS: 4 (0019, 106b) [Field of view in detector cells] SS: 852 (0019, 1070) [Value of back projection button] SS: 1 (0019, 1071) [Set if fatq estimates were used] SS: 0 (0019, 1072) [Z chan avg over views] DS: '0.0' (0019, 1073) [Avg of left ref chans over views] DS: '0.0' (0019, 1074) [Max left chan over views] DS: '0.0' (0019, 1075) [Avg of right ref chans over views] DS: '0.0' (0019, 1076) [Max right chan over views] DS: '0.0' (0019, 10da) [Reference channel used] SS: 0 (0019, 10db) [Back projector coefficient] DS: '0.0' (0019, 10dc) [Primary speed correction used] SS: 1 (0019, 10dd) [Overrange correction used] SS: 1 (0019, 10de) [Dynamic Z alpha value] DS: '0.0' (0020, 000d) Study Instance UID UI: 1.3.6.1.4.1.5962.1.2.1.20040119072730.12322 (0020, 000e) Series Instance UID UI: 1.3.6.1.4.1.5962.1.3.1.1.20040119072730.12322 (0020, 0010) Study ID SH: '1CT1' (0020, 0011) Series Number IS: '1' (0020, 0012) Acquisition Number IS: '2' (0020, 0013) Instance Number IS: '1' (0020, 0032) Image Position (Patient) DS: [-158.135803, -179.035797, -75.699997] (0020, 0037) Image Orientation (Patient) DS: [1.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000] (0020, 0052) Frame of Reference UID UI: 1.3.6.1.4.1.5962.1.4.1.1.20040119072730.12322 (0020, 0060) Laterality CS: '' (0020, 1040) Position Reference Indicator LO: 'SN' (0020, 1041) Slice Location DS: '-77.2040634155' (0020, 4000) Image Comments LT: 'Uncompressed' (0021, 0010) Private Creator LO: 'GEMS_RELA_01' (0021, 1003) [Series from which Prescribed] SS: 2 (0021, 1005) [Genesis Version - now] SH: '05' (0021, 1007) [Series Record checksum] UL: 1605775145 (0021, 1015) [Unknown] US: 24078 (0021, 1016) [Unknown] SS: 2 (0021, 1018) [Genesis version - Now] SH: '05' (0021, 1019) [Acq recon record checksum] UL: 750675506 (0021, 1037) [Screen Format] SS: 16 (0021, 104a) [Anatomical reference for scout] LO: '' (0021, 1090) [Tube focal spot position] SS: 7400 (0021, 1091) [Biopsy position] SS: 0 (0021, 1092) [Biopsy T location] FL: 0.0 (0021, 1093) [Biopsy ref location] FL: 0.0 (0023, 0010) Private Creator LO: 'GEMS_STDY_01' (0023, 1070) [Start time(secs) in first axial] FD: 862399761.111079 (0023, 1074) [No. of updates to header] SL: 1 (0023, 107d) [Indicates study has complete info SS: 0 (0025, 0010) Private Creator LO: 'GEMS_SERS_01' (0025, 1006) [Last pulse sequence used] SS: 0 (0025, 1007) [Images in Series] SL: 44 (0025, 1010) [Landmark Counter] SL: 0 (0025, 1011) [Number of Acquisitions] SS: 0 (0025, 1017) [Series Complete Flag] SL: 0 (0025, 1018) [Number of images archived] SL: 0 (0025, 1019) [Last image number used] SL: 4 (0025, 101a) [Primary Receiver Suite and Host] SH: '' (0027, 0010) Private Creator LO: 'GEMS_IMAG_01' (0027, 1006) [Image archive flag] SL: 1 (0027, 1010) [Scout Type] SS: 0 (0027, 101c) [Vma mamp] SL: 150 (0027, 101d) [Vma phase] SS: 1 (0027, 101e) [Vma mod] SL: 24 (0027, 101f) [Vma clip] SL: 129 (0027, 1020) [Smart scan ON/OFF flag] SS: 1 (0027, 1030) [Foreign Image Revision] SH: '' (0027, 1035) [Plane Type] SS: 2 (0027, 1040) [RAS letter of image location] SH: 'I' (0027, 1041) [Image location] FL: -77.20406341552734 (0027, 1042) [Center R coord of plane image] FL: -11.199999809265137 (0027, 1043) [Center A coord of plane image] FL: 9.699999809265137 (0027, 1044) [Center S coord of plane image] FL: -75.69999694824219 (0027, 1045) [Normal R coord] FL: 0.0 (0027, 1046) [Normal A coord] FL: 0.0 (0027, 1047) [Normal S coord] FL: -1.0 (0027, 1048) [R Coord of Top Right Corner] FL: -180.53579711914062 (0027, 1049) [A Coord of Top Right Corner] FL: 179.03579711914062 (0027, 104a) [S Coord of Top Right Corner] FL: -75.69999694824219 (0027, 104b) [R Coord of Bottom Right Corner] FL: -180.53579711914062 (0027, 104c) [A Coord of Bottom Right Corner] FL: -159.63580322265625 (0027, 104d) [S Coord of Bottom Right Corner] FL: -75.69999694824219 (0027, 1050) [Scan Start Location] FL: -63.19999694824219 (0027, 1051) [Scan End Location] FL: -116.20304870605469 (0027, 1052) [RAS letter for side of image] SH: 'L' (0027, 1053) [RAS letter for anterior/posterior] SH: 'A' (0027, 1054) [RAS letter for scout start loc] SH: 'I' (0027, 1055) [RAS letter for scout end loc] SH: 'I' (0028, 0002) Samples per Pixel US: 1 (0028, 0004) Photometric Interpretation CS: 'MONOCHROME2' (0028, 0010) Rows US: 128 (0028, 0011) Columns US: 128 (0028, 0030) Pixel Spacing DS: [0.661468, 0.661468] (0028, 0100) Bits Allocated US: 16 (0028, 0101) Bits Stored US: 16 (0028, 0102) High Bit US: 15 (0028, 0103) Pixel Representation US: 1 (0028, 0120) Pixel Padding Value SS: -2000 (0028, 1052) Rescale Intercept DS: '-1024.0' (0028, 1053) Rescale Slope DS: '1.0' (0029, 0010) Private Creator LO: 'GEMS_IMPS_01' (0029, 1004) [Lower range of Pixels1] SL: 0 (0029, 1005) [Lower range of Pixels1] DS: '0.0' (0029, 1006) [Lower range of Pixels1] DS: '0.0' (0029, 1007) [Lower range of Pixels1] SL: 87 (0029, 1008) [Lower range of Pixels1] SH: '' (0029, 1009) [Lower range of Pixels1] SH: '' (0029, 100a) [Lower range of Pixels1] SS: 764 (0029, 1026) [Version of the hdr struct] SS: 2 (0029, 1034) [Advantage comp. Overflow] SL: 0 (0029, 1035) [Advantage comp. Underflow] SL: 0 (0043, 0010) Private Creator LO: 'GEMS_PARM_01' (0043, 1010) [Window value] US: 400 (0043, 1011) [Total input views] US: 10431 (0043, 1012) [X-ray chain] SS: [14, 2, 3] (0043, 1013) [Decon kernel parameters] SS: [107, 21, 4, 2, 20] (0043, 1014) [Calibration parameters] SS: [4, 4, 5] (0043, 1015) [Total output views] SS: 10431 (0043, 1016) [Number of overranges] SS: 0 (0043, 1017) [IBH image scale factors] DS: '0.095' (0043, 1018) [BBH coefficients] DS: [0.085000, 1.102000, 0.095000] (0043, 1019) [Number of BBH chains to blend] SS: 350 (0043, 101a) [Starting channel number] SL: 7 (0043, 101b) [Ppscan parameters] SS: [0, 0, 0, 0, 0] (0043, 101c) [GE image integrity] SS: 0 (0043, 101d) [Level value] SS: 40 (0043, 101e) [Delta Start Time [msec]] DS: '2.0' (0043, 101f) [Max overranges in a view] SL: 0 (0043, 1020) [Avg overranges all views] DS: '0.0' (0043, 1021) [Corrected after glow terms] SS: 0 (0043, 1025) [Reference channels] SS: [1, 2, 3, 748, 749, 750] (0043, 1026) [No views ref chans blocked] US: [0, 1, 1, 0, 0, 0] (0043, 1027) [Scan Pitch Ratio] SH: '/1.0:1' (0043, 1028) [Unique image iden] OB: Array of 80 elements (0043, 1029) [Histogram tables] OB: Array of 2068 elements (0043, 102a) [User defined data] OB: Array of 40 elements (0043, 102b) [Private Scan Options] SS: [4, 4, 0, 0] (0043, 1031) [Recon Center Coordinates] DS: [-11.200000, 9.700000] (0043, 1040) [Trigger on position] FL: 178.07992553710938 (0043, 1041) [Degree of rotation] FL: 3816.219482421875 (0043, 1042) [DAS trigger source] SL: 0 (0043, 1043) [DAS fpa gain] SL: 0 (0043, 1044) [DAS output source] SL: 1 (0043, 1045) [DAS ad input] SL: 1 (0043, 1046) [DAS cal mode] SL: 3 (0043, 1047) [DAS cal frequency] SL: -1 (0043, 1048) [DAS reg xm] SL: 1 (0043, 1049) [DAS auto zero] SL: 0 (0043, 104a) [Starting channel of view] SS: 1 (0043, 104b) [DAS xm pattern] SL: 0 (0043, 104c) [TGGC trigger mode] SS: 0 (0043, 104d) [Start scan to X-ray on delay] FL: 0.0 (0043, 104e) [Duration of X-ray on] FL: 10.60060977935791 (7fe0, 0010) Pixel Data OW: Array of 32768 elements (fffc, fffc) Data Set Trailing Padding OB: Array of 126 elements
Matplotlib's imshow() function </b> cmap='gray' means gray color mappings for each value
import matplotlib.pyplot as plt
# Draw the image in grayscale
plt.imshow(im_ct, cmap='gray')
# Render the image
plt.show()
ImageIO's volread() function
Load multi-dimensional datasets and create 3D volumes from a folder of images
# Path to the folder containing the DICOM files
folder_path = 'Data/Brain/SE000001/MR000000'
# Load the brain data folder
vol = imageio.volread(folder_path)
Reading DICOM (examining files): 1/27 files (3.7%27/27 files (100.0%) Found 1 correct series. Reading DICOM (loading data): 27/27 (100.0%)
Slicing 3D and 4D images into many 2D frames
# Plot the first slice of the volume (2D frame)
plt.imshow(vol[0, :, :], cmap='gray')
# Render the image
plt.show()
Slicing along different axes
# Compute aspect ratios
d0, d1, d2 = vol.meta['sampling']
asp1 = d0 / d2
# Plot the other view slice of the volume
plt.imshow(vol[:, 128, :], cmap='gray',aspect=asp1)
# Render the image
plt.show()
Shifting or moving an image from one place to another within a picture.
import scipy.ndimage as ndi
#formatting method for the plots in this notebook:
def format_and_render_plot():
'''
Custom function to simplify common formatting operations for exercises. Operations include:
1. Turning off axis grids.
2. Calling `plt.tight_layout` to improve subplot spacing.
3. Calling `plt.show()` to render plot.
'''
fig = plt.gcf()
for ax in fig.axes:
ax.axis('off')
plt.tight_layout()
plt.show()
# save the middle slice as separat image
middle_slice = vol.shape[0] // 2 # // is floor division
im = vol[middle_slice,:,:]
# Translate the brain towards the center
xfm = ndi.shift(im, shift=(10, 15))
# Plot the original and adjusted images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('translated', fontweight ="bold")
format_and_render_plot()
Rotating an image from its center by the specified degrees from the right horizontal axis.
# Rotate the shifted image
xfm = ndi.rotate(xfm, angle=-30, reshape=False)
# Plot the original and transformed images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('rotated', fontweight ="bold")
format_and_render_plot()
Scaling (changing size), rotation, translation (shifting), and shearing (changing angles) using matrix multiplication
# Define the affine transform matrix
mat = np.array([[0.8, -0.4, 90], [0.4, 0.8, -6.0], [0, 0, 1]])
print(mat)
# Apply the affine transform matrix to image data
xfm = ndi.affine_transform(im, mat)
# Plot the original and transformed images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('transformed', fontweight ="bold")
format_and_render_plot()
[[ 0.8 -0.4 90. ] [ 0.4 0.8 -6. ] [ 0. 0. 1. ]]
# Resample image
im_dn = ndi.zoom(im, zoom=0.25)
im_up = ndi.zoom(im, zoom=4.00)
# Plot the images
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im_dn, cmap='gray')
axes[0].set_title('downsampled', fontweight ="bold")
axes[1].imshow(im_up, cmap='gray')
axes[1].set_title('upsampled', fontweight ="bold")
format_and_render_plot()
Estimating new pixel intensities after an image transformation
# Upsample "im" by a factor of 4
up0 = ndi.zoom(im, zoom=512/128, order=0)
up5 = ndi.zoom(im, zoom=512/128, order=5)
# Print original and new shape
print('Original shape:', im.shape)
print('Upsampled shape:', up5.shape)
# Plot close-ups of the new images
fig, axes = plt.subplots(1, 2)
axes[0].imshow(up0[128:256, 128:256], cmap='gray')
axes[0].set_title('no interpolation', fontweight ="bold")
axes[1].imshow(up5[128:256, 128:256], cmap='gray')
axes[1].set_title('with interpolation', fontweight ="bold")
format_and_render_plot()
Original shape: (256, 256) Upsampled shape: (1024, 1024)
Mean intensity differences between two images, with higher values indicating greater divergence
# Load the dcm file (image)
im1 = im
# Apply the affine transform matrix to image data
xfm = ndi.shift(im, shift=(-12, -8))
im2 = xfm
# Calculate image difference
err = im1 - im2
# Plot the difference
fig, axes = plt.subplots(1, 3)
axes[0].imshow(im1, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im2, cmap='gray')
axes[1].set_title('shifted', fontweight ="bold")
axes[2].imshow(err, cmap='seismic', vmin=-200, vmax=200)
axes[2].set_title('error map', fontweight ="bold")
format_and_render_plot()
# Calculate absolute image difference
abs_err = np.absolute(im1 - im2)
# Plot the difference
fig, axes = plt.subplots(1, 3)
axes[0].imshow(im1, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im2, cmap='gray')
axes[1].set_title('shifted', fontweight ="bold")
axes[2].imshow(abs_err, cmap='seismic', vmin=-200, vmax=200)
axes[2].set_title('error map (absolute)', fontweight ="bold")
format_and_render_plot()
# Calculate mean absolute error
mean_abs_err = np.mean(np.abs(im1 - im2))
print('MAE:', mean_abs_err)
MAE: 25730.311309814453
Structural similarity between two images perceived changes in:
structural information, luminance, and contrast that humans typically notice.
from skimage.metrics import structural_similarity as ssim
# Compute SSIM between im and xfm
data_range=xfm.max() - xfm.min()
ssim_index = ssim(im, xfm, data_range=data_range)
print(f'SSIM: {ssim_index}')
SSIM: 0.47951331548119075
The number of pixels filled in both images (the intersection) out of the number of pixels filled in either image (the union).
def intersection_of_union(im1, im2):
i = np.logical_and(im1, im2)
u = np.logical_or(im1, im2)
return i.sum() / u.sum()
# Try some other paramters by yourself
xfm = ndi.shift(im, shift=(-10, -10))
xfm = ndi.rotate(xfm, angle=-15, reshape=False)
intersection_of_union(xfm, im)
0.8135092272202998
The number of pixels filled in both images (the intersection) out of the number of pixels filled in either image (the union).
import SimpleITK as sitk
import Utilities.gui
img1 = sitk.ReadImage('../ImageAnalysis/Data/Brain/training_001_mr_T1.mha') #MRI image
img2_original = sitk.ReadImage('../ImageAnalysis/Data/Brain/training_001_ct.mha') #CT image
img2 = sitk.Resample(img2_original, img1)
# Obtain foreground masks for the two images using Otsu thresholding, we use these later on.
msk1 = sitk.OtsuThreshold(img1,0,1)
msk2 = sitk.OtsuThreshold(img2,0,1)
Utilities.gui.MultiImageDisplay(image_list = [img1, img2],
title_list = ['image1', 'image2'],
figure_size=(9,3))
<Utilities.gui.MultiImageDisplay at 0x2a5c2f58280>
# now we load test data files (CT and MR)
ct_file = get_testdata_file("CT_small.dcm")
# Load the CT file (image)
im = imageio.imread(ct_file)
C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/2938429760.py:5: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly. im = imageio.imread(ct_file)
# min-max normalisation
im_old = im # save original image
print('Min. value before normalization:', im_old.min())
print('Max value before normalization:', im_old.max())
im = (im - im.min()) / (im.max() - im.min()) # normalise to 0-1 range
print('Min. value after normalization:', im.min())
print('Max value after normalization:', im.max())
Min. value before normalization: -896 Max value before normalization: 1167 Min. value after normalization: 0.0 Max value after normalization: 1.0
def format_and_render_plot(axis=False, legend=False):
'''
Custom function to simplify common formatting operations for exercises. Operations include:
1. Turning off axis grids and legends, if not explicitly requested.
2. Calling `plt.tight_layout` to improve subplot spacing.
3. Calling `plt.show()` to render plot.
'''
fig = plt.gcf()
for ax in fig.axes:
if not axis:
ax.axis('off')
if legend:
ax.legend(loc='center right')
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].imshow(im_old, cmap='gray')
axes[0].set_title('without normalization', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].set_title('with normalization', fontweight ="bold")
format_and_render_plot()
# Create a 256-bin histogram, binned at each possible value
hist = ndi.histogram(im, min=im.min(), max=im.max(), bins=256)
# Create a cumulative distribution function
cdf = hist.cumsum() / hist.sum()
# Plot the histogram and CDF
fig, axes = plt.subplots(2, 1, sharex=True) # sharex=True shares the x-axis between the top and bottom subplot
axes[0].plot(hist, label='Histogram') # label='Histogram' labels this line as "Histogram" for the legend
axes[1].plot(cdf, label='CDF') # label='CDF' labels this line as "CDF" for the legend
format_and_render_plot(axis=True, legend=True) # axis=True turns on axis grids for the plot; legend=True turns on the legend
Binary arrays for removing or selecting specific parts of an image by applying one or more logical operations to an image.
# Create soft tissue (st) and bone masks
# Try different intervals by yourself
mask_st = (im >= 0.2) & (im < 0.52) # soft tissue mask
mask_bone = im >= 0.53 # bone mask
# Plot the skin (0) and bone (1) masks
fig, axes = plt.subplots(1,2)
axes[0].imshow(mask_st, cmap='gray')
axes[0].set_title('soft tissue map', fontweight ="bold")
axes[1].imshow(mask_bone, cmap='gray')
axes[1].set_title('bone map', fontweight ="bold")
format_and_render_plot()
Although masks are binary, they can be applied to images to filter out pixels where the mask is False.
NumPy's where() function is a flexible way of applying masks. It takes three arguments:
np.where(condition, x, y)
condition, x and y can be either arrays or single values. This allows you to pass through original image values while setting masked values to 0.
# Screen out non-bone pixels from "im"
im_bone = np.where(mask_bone, im, 0) #extract pixels from im which locate in mask_bone!
# Plot masked image and histogram
plt.imshow(im_bone, cmap='gray', label='Bone pixels')
<matplotlib.image.AxesImage at 0x2a5c6fd0520>
Imperfect masks can be tuned through the addition and subtraction of pixels. example methods for accomplishing these ends:
binary_dilation: Add pixels along edgesbinary_closing: Dilate then erode, "filling in" holes# Create and tune bone mask
mask_dilate = ndi.binary_dilation(mask_bone, iterations=5)
mask_closed = ndi.binary_closing(mask_bone, iterations=5)
# Plot masked images
fig, axes = plt.subplots(1,3)
axes[0].imshow(mask_bone)
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(mask_dilate)
axes[1].set_title('dilated', fontweight ="bold")
axes[2].imshow(mask_closed)
axes[2].set_title('closed', fontweight ="bold")
format_and_render_plot()
Transforming images based on intensity values surrounding a pixel
# Set filter weights
weights = [[0.11, 0.11, 0.11],
[0.11, 0.11, 0.11],
[0.11, 0.11, 0.11]]
# Convolve the image with the filter
im_filt = ndi.convolve(im, weights)
# Plot the images
fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_filt, cmap='gray')
axes[1].set_title('with filter', fontweight ="bold")
format_and_render_plot()
Improve the signal-to-noise ratio of an image by blurring out small variations in intensity. The Gaussian filter is excellent for this: it is a circular (or spherical) smoothing kernel that weights nearby pixels higher than distant ones.
# Smooth "im" with Gaussian filters
# Try different sigmas by yourself
im_s1 = ndi.gaussian_filter(im, sigma=1)
im_s3 = ndi.gaussian_filter(im, sigma=3)
# Draw bone masks of each image
fig, axes = plt.subplots(1,3)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_s1, cmap='gray')
axes[1].set_title('Gaussian (sigma=1)', fontweight ="bold")
axes[2].imshow(im_s3, cmap='gray')
axes[2].set_title('Gaussian (sigma=3)', fontweight ="bold")
format_and_render_plot()
Applying filters which pattern is a change in intensity along a plane.
# Set weights to detect vertical edges
weights = [[1, 0, -1],
[1, 0, -1],
[1, 0, -1]]
# Convolve "im" with filter weights
edges = ndi.convolve(im, weights)
# Draw the image in color
fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(edges, cmap='seismic', vmin=-0.8, vmax=0.8)
axes[1].set_title('edge filter', fontweight ="bold")
format_and_render_plot()
# Load the directory (volume)
folder_path = 'Data/Brain/SE000001/MR000000'
# Load the volume
vol = imageio.volread(folder_path)
print(vol.shape)
# save the middle slice as separat image
middle_slice = vol.shape[0] // 2 # // is floor division
im = vol[middle_slice,:,:]
# min-max normalisation
im_old = im # save original image for later
im = (im - im.min()) / (im.max() - im.min()) # normalise the image, range 0 - 1
Reading DICOM (examining files): 1/27 files (3.7%27/27 files (100.0%) Found 1 correct series. Reading DICOM (loading data): 27/27 (100.0%) (27, 256, 256)
Apply a mask and fill any gaps in the mask
# Smooth intensity values
im_filt = ndi.median_filter(im, size=3) # size = 3 means 3x3x3 neighbourhood
# Select high-intensity pixels
mask_start = np.where(im_filt > 0.3, 1, 0) # mask_start is a boolean array
mask = ndi.binary_closing(mask_start) # fill holes
# Label the objects in "mask"
labels, nlabels = ndi.label(mask) # labels: each object has a unique number, nlabels: number of objects
# Create a `labels` overlay
overlay = np.where(labels > 0, labels, np.nan)
# Use imshow to plot the overlay
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray') # show image first
axes[1].imshow(overlay, cmap='rainbow', alpha=0.6) # show overlay second; alpha controls transparency
axes[1].set_title('with segmentation', fontweight ="bold")
format_and_render_plot()
Pick up whole sets of pixels at a time
# Label the image "mask"
labels, nlabels = ndi.label(mask)
# Select brain label
brain_val = 2 # 2 is the label for the brain (see plot above)
brain_mask = np.where(labels == brain_val, 1, np.nan) # create brain mask
# Overlay selected label
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].imshow(brain_mask, cmap='rainbow', alpha=0.6)
axes[1].set_title('with segmentation', fontweight ="bold")
format_and_render_plot()
Crop images so that they only include the object of interest.
brain_mask=brain_mask.astype(np.int64)
# Find bounding box
bboxes =ndi.find_objects(brain_mask) # returns a list of tuples, each tuple has 3 slices
print('Number of objects:', len(bboxes)) # number of objects found
print('Indices for first box:', bboxes[0]) # print indices for first box
# Crop to index 0
im_brain = im[bboxes[0]] # crop the original image to the bounding box
# Plot the cropped image
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_brain, cmap='gray')
axes[1].set_title('cropped', fontweight ="bold")
format_and_render_plot()
Number of objects: 1 Indices for first box: (slice(50, 234, None), slice(61, 193, None))
Tailor measurements to specific sets of pixels:
# Variance for all pixels
var_all = ndi.variance(vol, labels=None, index=None) # labels=None means all pixels, index=None means all objects
print('All pixels:', var_all)
# Variance for labeled pixels
var_labels = ndi.variance(vol, labels, index=None) # labels=labels means only labeled pixels, index=None means all objects
print('Labeled pixels:', var_labels)
# Variance for each object
var_objects = ndi.variance(vol, labels, index=[1,2]) # labels=labels means only labeled pixels, index=[1,2] means only objects 1 and 2
print('Brain matter:', var_objects[0])
print('Other tissue:', var_objects[1])
All pixels: 15281.189220274951 Labeled pixels: 18188.560688187546 Brain matter: 13466.402775068886 Other tissue: 16319.82340391041
Wide distribution of intensity values and more variance of multiple tissue types
# Create histograms for selected pixels
hist1 = ndi.histogram(vol, min=0, max=255, bins=256)
hist2 = ndi.histogram(vol, 0, 255, 256, labels=labels)
hist3 = ndi.histogram(vol, 0, 255, 256, labels=labels, index=1)
# Plot the histogram and CDF
fig, axes = plt.subplots(3, 1, sharex=True) # sharex=True shares the x-axis between the top and bottom subplot
axes[0].plot(hist1 / hist1.sum(), label='All pixels')
axes[1].plot(hist2 / hist2.sum(), label='All labeled pixels')
axes[2].plot(hist3 / hist3.sum(), label='Brain matter')
format_and_render_plot(axis=True, legend=True) # axis=True turns on axis grids for the plot; legend=True turns on the legend
the distance from each pixel to a given point
# Calculate distances
brain = np.where(labels == 2, 1, 0) # create brain mask
dists = ndi.distance_transform_edt(brain, sampling=vol.meta['sampling'][1:3]) # calculate distances
# Report on distances
print('Max distance (mm):', ndi.maximum(dists))
print('Max location:', ndi.maximum_position(dists))
# Plot overlay of distances
overlay = np.where(dists > 0, dists, np.nan) # create overlay
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].imshow(overlay, cmap='hot', alpha=0.75)
axes[1].set_title('with distances', fontweight ="bold")
format_and_render_plot()
Max distance (mm): 12.1875 Max location: (133, 175)
The coordinates for the point of an object with higher values pulling the center closer to it.
# Extract centers of mass for objects 1 and 2
coms = ndi.center_of_mass(vol, labels, index=[1,2])
print('Label 1 center:', coms[0])
print('Label 2 center:', coms[1])
# Add marks to plot
for c0, c1, c2 in coms:
plt.scatter(c2, c1, s=100, marker='o')
plt.show()
Label 1 center: (8.856682792208547, 49.44438084201855, 125.03554483621987) Label 2 center: (11.713445688803832, 119.3292456288401, 129.36906622106693)